#pragma once
#include "hdev.h"


#include <acml.h>
#include <math.h>
#include <complex>

#ifdef __CUDACC__
#  include <cublas.h>
#endif

template <typename float_t> struct mathfns;

template <> struct mathfns<float>{
	static inline _hdev_ float sqrt (float x){
		return ::sqrtf(x);
	}
};

template <> struct mathfns<double>{
	static inline _hdev_ double sqrt (double x){
		return ::sqrt(x);
	}
};

template <typename float_t>
struct cplx_t
{

	float_t real, imag;

	inline _hdev_ cplx_t (){}
	inline _hdev_ cplx_t (float_t r, float_t i) : real(r), imag(i) {}

	inline cplx_t (complex const& x) : real(x.real), imag(x.imag) {}

	#ifdef __CUDACC__
	inline _hdev_ cplx_t (cuComplex const& x) : real(x.x), imag(x.y) {}
	#endif

	template <typename float_2>
	inline cplx_t (std::complex<float_2> const& x){
		this->real = std::real(x);
		this->imag = std::imag(x);
	}

	inline operator complex() const {
		return compose_complex (real, imag);
	}

	template <typename float_2>
	inline operator std::complex<float_2>() const {
		return std::complex<float_2>(real, imag);
	}

	inline operator doublecomplex() const {
		return compose_doublecopmlex (real, imag);
	}

	inline _hdev_ cplx_t conj() const {
		return cplx_t (real, -imag);
	}

	inline _hdev_ float_t abs() const {
		return mathfns<float_t>::sqrt (real*real + imag*imag);
	}

	inline _hdev_ float_t abs2() const {
		return real*real + imag*imag;
	}

	inline _hdev_ cplx_t operator+ (cplx_t const& y) const {
		return cplx_t (real + y.real, imag + y.imag);		
	}
	inline _hdev_ cplx_t operator- (cplx_t const& y) const {
		return cplx_t (real - y.real, imag - y.imag);		
	}

	inline _hdev_ cplx_t operator- () const {
		return cplx_t (-real, -imag);
	}

	inline _hdev_ cplx_t operator* (cplx_t const& y) const {
		return cplx_t (real*y.real - imag*y.imag, real*y.imag + imag*y.real);
	}

	inline _hdev_ cplx_t operator* (float alpha) const {
		return cplx_t (alpha * real, alpha * imag);
	}

	inline _hdev_ cplx_t operator* (double alpha) const {
		return cplx_t (alpha * real, alpha * imag);
	}

	inline _hdev_ cplx_t& operator= (cplx_t const& y) {
		real = y.real; imag = y.imag;
		return *this;
	}

	inline _hdev_ cplx_t& operator= (float y) {
		real = y; imag = 0;
		return *this;
	}

	inline _hdev_ cplx_t& operator= (double y) {
		real = y; imag = 0;
		return *this;
	}

	inline _hdev_ cplx_t& operator+= (cplx_t const& y){
		real += y.real; imag += y.imag;
		return *this;
	}

	inline _hdev_ cplx_t& operator-= (cplx_t const& y){
		real -= y.real; imag -= y.imag;
		return *this;
	}

	inline _hdev_ cplx_t& operator*= (cplx_t const& y) {
		float_t r = real*y.real - imag*y.imag;
		float_t i = real*y.imag + imag*y.real;
		real = r; imag = i;
		return *this;
	}

	inline _hdev_ cplx_t& operator*= (float alpha) {
		real *= alpha; imag *= alpha;
		return *this;
	}

	inline _hdev_ cplx_t& operator*= (double alpha) {
		real *= alpha; imag *= alpha;
		return *this;
	}

};

template <typename float_t>
static inline _hdev_ cplx_t<float_t> conj (cplx_t<float_t> const& x){
	return x.conj();
}

template <typename float_t>
static inline _hdev_ float_t abs (cplx_t<float_t> const& x){
	return x.abs();
}

template <typename float_t>
static inline _hdev_ float_t abs2 (cplx_t<float_t> const& x){
	return x.abs2();
}

template <typename float_t>
static inline _hdev_ cplx_t<float_t> operator* (float alpha,
                                         cplx_t<float_t> const& x){
	return x*alpha;
}

template <typename float_t>
static inline _hdev_ cplx_t<float_t> operator* (double alpha,
                                         cplx_t<float_t> const& x){
	return x*alpha;
}


typedef cplx_t<float> cplxf;
typedef cplx_t<double> cplxd;
